##############################
# Replication code for Harbers & Ingram, "Geo-Nested Analysis" (Political Analysis)
# Data from Messner et al. 2000 and Baller et al. 2001 (see README file and citations in paper)
# Underlying data used in illustration of method are publicly avaiable at: http://spatial.uchicago.edu/sample-data, and
# are also included in replication materials to facilitate re-analysis
# Last updated: 2016-12-28 (MI)
##############################

# note for new R users: pound sign ('#') identifies a comment line that is not run
# code below is extensively commented for ease of replication by both new and experienced R users

#################################
# set working directory (can also do this manually from File pulldown menu)
# NOTE: this is set to my own working directory; select your own
#################################

# if don't know working directory, use getwd()
getwd()

# if running on cluster or network computer:
path <- '/network/rit/lab/ingramlab/geonested'  # again, this is my working directory; set your own
# or, if running from local computer (not secure shell)
# path <- setwd('A:/geonested')    # where 'A' is mapped network drive
# or, if using directory on disk drive
# path <- 'C:/Users/MI122167/Dropbox/Met Matt/Mixed Methods Paper/PA Submission/replication'
setwd(path)
# path <- 'C:/Users/mi122167/Dropbox/SUNY Albany/Spatial Analysis/Data and Replications/HarbersIngram2016-APSA/RPOS 619 GWR-SL'
dir()
# can also check directory contents with: list.files()

############################
# install necessary packages
############################
# install.packages()
#install.packages("spdep")
#install.packages("graphicsQC")
#install.packages("prettyR")
#install.packages("RColorBrewer")
#install.packages("GWmodel")
#install.packages("car")
#install.packages("spgwr")
#install.packages("gridExtra")
#install.packages("sp")
#install.packages("classInt")
#install.packages("rgeos")
#install.packages("maptools")
#install.packages("spdep")
#install.packages("RgoogleMaps")
#install.packages("rgdal")
#install.packages("weights")
#install.packages("INLA")
#install.packages("xts")
#install.packages("GeoXp")
#install.packages("spacetime")
#install.packages("plm")
#install.packages("ggplot2")
#install.packages("INLA")
#install.packages("maps")
#install.packages("plm")
#install.packages("INLA")
#install.packages("latticeExtra")
#install.packages("plyr")
#install.packages("graphicsQC")
#install.packages("parallel")
#install.packages("scales")
#install.packages("cowplot")
#install.packages("sphet")
#install.packages("texreg")

########################
# note for new R users: installing is not enough; need to open library with library() command
# note: use 'library' command instead of 'require'; 'require' ends up calling 'library' command, but does not issue error warnings
########################

# open libraries 
#library(sp)
library(spdep)  # loads sp, along with other dependencies
#library(classInt)
#library(rgeos)
library(maptools)
gpclibPermit()   # if get error msg from maptools, use this to enable gpclib for maptools
#library(RgoogleMaps)
#library(rgdal)
library(ggplot2)
library(weights)
library(car)
#library(GeoXp)
#library(spacetime)
#library(xts)
#library(INLA)
#library(maps)
#library(plm)
library(RColorBrewer)
library(latticeExtra)
library(prettyR)
library(gridExtra)
library(plyr)
library(GWmodel)
#library(graphicsQC)
library(spgwr)
#library(parallel)
library(scales)
#library(cowplot)
#library(sphet)
library(texreg)

# check to make sure all packages and dependencies loaded correctly
sessionInfo()

##################
# READING/IMPORTING AND WRITING/EXPORTING DATA (shapefile, .shp)
##################

# get shapefile info
# getinfo.shape
# getinfo.shape("NAT.SHP")
# load shapefile and attaching data
dir()
setwd(paste0(path, '/shapefiles-data', sep=""))

shp <- readShapeSpatial("NAT.SHP", ID="FIPSNO")
# note: sometimes dbf file will not load if extension is .DBF; if this happens, try renaming to lower case .dbf

shp.so <- readShapeSpatial("nat_south2.SHP", ID="FIPSNO")
shpns <- readShapeSpatial("nat_nonsouth2.SHP", ID="FIPSNO")

shp.st <- readShapeSpatial("tl_2014_us_state_contig.shp")
shp.st.so <- readShapeSpatial("tl_2014_us_state_contig_south.shp")
shp.st.ns <- readShapeSpatial("tl_2014_us_state_contig_nonsouth.shp")

#plot(shp)
#plot(shp.st, border="red", add=TRUE)

#summary(shp)
data <- data.frame(shp)
data.so<-data.frame(shp.so)
data.ns <- data.frame(shpns)
#names(data)
#head(data)
#detach()
#attach(data.so)

#################
# correlations table, if desired
#################

#desc.vars <- cbind(FIPSNO, HR90, RD90, PS90, MA90, DV90, UE90, SOUTH)
#describe(desc.vars, num.desc = c("valid.n", "min", "mean", "max", "sd")) # package prettyR
#cor(desc.vars[,2:8])   #correlation matrix excluding the ID variable


##############
# We need a grid reference for the location of each observation.
# we create this location using the coordinates of the centroids:

# need to merge X-Y coords with shapefile and then re-load shapefile and attach data
# XYgridtable <- cbind(XCNTLONG, YCNTLAT)
# or simply:
coords <- coordinates(shp)
coords.so <- coordinates(shp.so)
coords.ns <- coordinates(shpns)

###############
#
# Spatial Weights
#
###############

# Here, we might only calculate 10 nearest neighbors (10nnb) since that is what is used by Baller et al
# However, we will calculate that W (w10), as well as first-order rook (r1) and queen (q1) contiguity matrices 
# row.names(shp)

# for full data set (all 3085 counties)
IDs<-row.names(as(shp, "data.frame"))
w_10nnb<-knn2nb(knearneigh(coords, k=10), row.names=IDs)

# make symmetric
wsym_10nnb <- make.sym.nb(w_10nnb)
wsym_10nnb

# and make this a list objects
listwsym_10nnb <- nb2listw(wsym_10nnb, style="W")
listwsym_10nnb

# create shortcut to symmetric list objects
w10 <- listwsym_10nnb
w10 


# for only south (1412 counties)
IDs<-row.names(as(shp.so, "data.frame"))
w_10nnb.so<-knn2nb(knearneigh(coords.so, k=10), row.names=IDs)

# make symmetric
wsym_10nnb.so <- make.sym.nb(w_10nnb.so)
wsym_10nnb.so

# second-order 10nnb w
wsym_10nnb2.so <- nblag(wsym_10nnb.so, maxlag=2)
wsym_10nnb2.so
# not that this is a list object containing two nb objects
wsym_10nnb2.so[[1]]  # this is first order w, equivalent to wsym_10nnb
wsym_10nnb2.so[[2]]  # this is second order w

# if you wanted cumulative lags for all orders within nblag object (here, wsym_10nnb2):
# w.cumul <- nblag_cumul(wsym_10nnb2)

# and make these list objects
listwsym_10nnb.so <- nb2listw(wsym_10nnb.so, style="W")
listwsym_10nnb.so
listwsym_10nnb2.so <- nb2listw(wsym_10nnb2.so[[2]], style="W")
listwsym_10nnb2.so

# create shortcut to symmetric list objects
w10s <- listwsym_10nnb.so
w10s 

w10.2s <- listwsym_10nnb2.so
w10.2s

# plot, if desired
#par(mfrow=c(2,1)) 
#plot(w10, coords, col="red")
#plot(w10.2, coords, col="blue")
#par(mfrow=c(1,1))

#plot(shp)
#plot(w10, coords, col="red", add=T)
#plot(w10.2, coords, col="blue", add=T)

# rook 1
w_nbr1<-poly2nb(shp.so, queen=FALSE)
# make symmetric
wsym_nbr1 <- make.sym.nb(w_nbr1)
wsym_nbr1

# second-order 10nnb w
wsym_nbr1.2 <- nblag(wsym_nbr1, maxlag=2)
wsym_nbr1.2

# and make these list objects
listwsym_nbr1 <- nb2listw(wsym_nbr1, style="W")
listwsym_nbr1
listwsym_nbr1.2 <- nb2listw(wsym_nbr1.2[[2]], style="W")
listwsym_nbr1.2
# create shortcuts
wr1s <- listwsym_nbr1
wr1.2s <- listwsym_nbr1.2

wr1s
wr1.2s

# plot
#par(mfrow=c(2,1)) 
#plot(wr1, coords, col="red")
#plot(wr1.2, coords, col="blue")
#par(mfrow=c(1,1))

# queen 1
w_nbq1<-poly2nb(shp.so, queen=TRUE)
# make symmetric
wsym_nbq1 <- make.sym.nb(w_nbq1)
wsym_nbq1

# second-order 10nnb w
wsym_nbq1.2 <- nblag(wsym_nbq1, maxlag=2)
wsym_nbq1.2

# and make these list objects
listwsym_nbq1 <- nb2listw(wsym_nbq1, style="W")
listwsym_nbq1
listwsym_nbq1.2 <- nb2listw(wsym_nbq1.2[[2]], style="W")
listwsym_nbq1.2
# create shortcuts
wq1s <- listwsym_nbq1
wq1.2s <- listwsym_nbq1.2

wq1s
wq1.2s

# plot
#par(mfrow=c(2,1)) 
#plot(wq1, coords, col="red")
#plot(wq1.2, coords, col="blue")
#par(mfrow=c(1,1))

setwd(path)
dir()
save.image("harbersingram20161223.RData")
#load("harbersingram20161116.RData")

# non-south

coords.ns <- coordinates(shpns)
IDs<-row.names(as(shpns, "data.frame"))
w_10nnb.ns <- knn2nb(knearneigh(coords.ns, k=10), row.names=IDs)

# make symmetric
wsym_10nnb.ns <- make.sym.nb(w_10nnb.ns)
wsym_10nnb.ns

# second-order 10nnb w
wsym_10nnb2.ns <- nblag(wsym_10nnb.ns, maxlag=2)
wsym_10nnb2.ns
# not that this is a list object containing two nb objects
wsym_10nnb2.ns[[1]]  # this is first order w, equivalent to wsym_10nnb
wsym_10nnb2.ns[[2]]  # this is second order w

# if you wanted cumulative lags for all orders within nblag object (here, wsym_10nnb2):
# w.cumul <- nblag_cumul(wsym_10nnb2)

# and make these list objects
listwsym_10nnb.ns <- nb2listw(wsym_10nnb.ns, style="W")
listwsym_10nnb.ns
listwsym_10nnb2.ns <- nb2listw(wsym_10nnb2.ns[[2]], style="W")
listwsym_10nnb2.ns

# create shortcut to symmetric list objects
w10ns <- listwsym_10nnb.ns
w10ns 

w10.2ns <- listwsym_10nnb2.ns
w10.2ns

# plot
#par(mfrow=c(2,1)) 
#plot(w10, coords, col="red")
#plot(w10.2, coords, col="blue")
#par(mfrow=c(1,1))

#plot(shp)
#plot(w10, coords, col="red", add=T)
#plot(w10.2, coords, col="blue", add=T)

# rook 1
w_nbr1<-poly2nb(shpns, queen=FALSE)
# make symmetric
wsym_nbr1 <- make.sym.nb(w_nbr1)
wsym_nbr1

# second-order 10nnb w
wsym_nbr1.2 <- nblag(wsym_nbr1, maxlag=2)
wsym_nbr1.2

# and make these list objects
listwsym_nbr1 <- nb2listw(wsym_nbr1, style="W")
listwsym_nbr1
listwsym_nbr1.2 <- nb2listw(wsym_nbr1.2[[2]], style="W")
listwsym_nbr1.2
# create shortcuts
wr1ns <- listwsym_nbr1
wr1.2ns <- listwsym_nbr1.2

wr1ns
wr1.2ns

# plot
#par(mfrow=c(2,1)) 
#plot(wr1, coords, col="red")
#plot(wr1.2, coords, col="blue")
#par(mfrow=c(1,1))

# queen 1
w_nbq1<-poly2nb(shpns, queen=TRUE)
# make symmetric
wsym_nbq1 <- make.sym.nb(w_nbq1)
wsym_nbq1

# second-order 10nnb w
wsym_nbq1.2 <- nblag(wsym_nbq1, maxlag=2)
wsym_nbq1.2

# and make these list objects
listwsym_nbq1 <- nb2listw(wsym_nbq1, style="W")
listwsym_nbq1
listwsym_nbq1.2 <- nb2listw(wsym_nbq1.2[[2]], style="W")
listwsym_nbq1.2
# create shortcuts
wq1ns <- listwsym_nbq1
wq1.2ns <- listwsym_nbq1.2

wq1ns
wq1.2ns

setwd(path)
dir()
save.image("harbersingram20161223.RData")

###################
#
# Table 1
#
###################

##################
# OLS global model, with regional dummy
##################
ols1 <- lm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + SOUTH, data = data)
#summary(ols1)
#AIC(ols1)

##################
# OLS global model, without regional dummy
##################
ols2 <- lm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90, data = data)
#summary(ols2)
#AIC(ols2)

# str(ols2)
# will need residual from ols2 later; ols2$residuals

##################
# SLM in south
##################

slm <- lagsarlm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90, data = data.so, listw = w10s, type="lag", method="eigen")

##################
# SEM in non-south
##################

sem <- errorsarlm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90, data = data.ns, listw = w10ns, etype="error", method="eigen")

###################
# Table 1
###################

# output results to file (csv or doc)
ols1sum <- summary(ols1)
ols2sum <- summary(ols2)
slmsum <- summary(slm)
semsum <- summary(sem)
msum <- c(ols1sum, ols2sum, slmsum, semsum)
setwd(paste0(path, "/output", sep=""))

# clean output of table 1 (html file, which can be copied and pasted into Word doc)
# note 1: htmlreg command output file with extension .doc; this table can be opened and edited in Word or text editor;
# can replace .doc with .htm to generate html file for web browser;
# note 2: to output table for latex, replace 'htmlreg' command in syntax below with 'texreg', and 
# in file output option replace file extension .htm with .tex

library(texreg)
htmlreg(list(ols1, ols2, slm, sem), 
          custom.model.names=c("1: OLS","2: OLS Without Regional Dummy", "3: Lag Model South", "4: Error Model Non-South"), 
          reorder.coef=c(2:9, 1),
          custom.coef.names = c("Constant", "Resource Dep. (rd90)", "Population Structure (ps90)", "Median Age (ma90)", 
                                "Divorce Rate (dv90)", "Unemployment Rate (ue90)", "South", "Rho", "Lambda"),
          digits=3,
          include.rmse = FALSE, include.nobs = TRUE, include.loglik = FALSE, include.lr = FALSE,
          custom.note = "Standard errors in parentheses. %stars.",
          file="table1.htm")

# to examine full, unstructured output
capture.output(msum, file="table1.csv")

setwd(path)

##########################################
#
# Figure 2: LISA values and cluster map
#
##########################################

# best practice: calculate 3 values (default with assumptions of equal variance, and two conservative estimates: saddlepoint and exact)
# if all close, then no problem using default; if not, then use one of more conservative methods
w<-w10

lisa <- localmoran(shp$HR90, w) # basic local moran with default settings with normalization assumption (default also means no probability value adjustments)

# for saddlepoint method, need nb object
nb <- w_10nnb
model <- lm(HR90 ~ 1, data=data) 
lisa_s <- localmoran.sad(model, nb=nb)  # if get error about different object lengths, make sure both model data and nb object are for full data (3085 counties)

# exact method
lisa_e <- localmoran.exact(model, nb=nb)

summary(lisa)
summary(lisa_s)
summary(lisa_e)

lisa
lisa_s
lisa_e

lisa_s.sum <- summary.localmoransad(lisa_s)
lisa_s.sum[1,]

# convert lisa objects to data frame objects
lisa <- as.data.frame(localmoran(shp$HR90, w))
lisa_s <- as.data.frame(localmoran.sad(lm(shp$HR90 ~ 1, shp), nb=nb, style="C")) # this works
lisa_e <- as.data.frame(localmoran.exact(lm(HR90 ~ 1, shp), nb=nb, style="C")) # this works

head(lisa)
head(lisa_s)
head(lisa_e)

# means should match global Moran's I
mean(lisa[,1])
lm.morantest(model, nb2listw(nb))
mean(lisa_s[,1]) # reports mean slightly higher than under normalization assumption
mean(lisa_e[,1]) # reports same mean as lisa, normalization assumption

# compare histograms of three LISA estimates
# prettier combination graph
hist1 <- hist(lisa[,"Pr(z > 0)"],breaks=20)
hist2 <- hist(lisa_s[,"Pr. (Sad)"],breaks=20)
hist3 <- hist(lisa_e[,"Pr. (exact)"],breaks=20)
par(mar=c(0,2,1,0)) # sets margin parameters for plot space
par(mfrow=c(3,1)) # sets graphing frame of 3 rows and 1 column
plot(hist1, border = "dark blue", col = "light blue",
     main = "LISA p values (norm)", xlab = "Pr")
plot(hist2, border = "dark blue", col = "light blue",
     main = "LISA p values (sad)", xlab = "Pr")
plot(hist3, border = "dark blue", col = "light blue",
     main = "LISA p values (exact)", xlab = "Pr")
par(mfrow=c(1,1)) # resets your graph space

# correlations matrix
desc.vars <- cbind(lisa[,"Ii"], lisa_s[,"Local Morans I"], lisa_e[,"Local Morans I"])
describe(desc.vars, num.desc = c("valid.n", "min", "mean", "max", "sd")) # package prettyR
cor(desc.vars)   #correlation matrix excluding the ID variable
# note: highly correlated

# correlations scatter plot
# simple version 1
pairs(~V1+V2+V3, data=lisa.data,
      main="Three LISA Values")

# car library (can also do others, e.g., lattice library)
library(car)
lisa.data <- as.data.frame(cbind(lisa[,1], lisa_s[,1], lisa_e[,1]))
scatterplotMatrix(~V1+V2+V3, data=lisa.data,
                  main="Three LISA Values")
# shows some divergence between (1) basic lisa and (2) both sad and exact at extreme negative values of LISA

# create LISA cluster identifiers, using saddlepoint data
DV <- shp$HR90
quadrant <- vector(mode="numeric",length=nrow(lisa_s))
cDV <- DV - mean(DV) 
lagDV <- lag.listw(w, DV)
clagDV <- lagDV - mean(lagDV)

plot(lagDV, DV)
plot(cDV, clagDV)
plot(ecdf(DV))
plot(ecdf(lagDV), add=TRUE, col.points="red", col.hor="red")


# LISA significance, using saddlepoint data
p <- lisa_s[,5]
quadrant <- vector(mode="numeric",length=nrow(lisa_s))
quadrant[cDV >0 & clagDV>0 & p<=.05] <- 1 
quadrant[cDV <0 & clagDV<0 & p<=.05] <- 2      
quadrant[cDV <0 & clagDV>0 & p<=.05] <- 3
quadrant[cDV >0 & clagDV<0 & p<=.05] <- 4
# non-significant obs will remain coded as zeroes (0)

# merge LISA of DV data into shapefile
# merge a single variable from table
shp$lisa_s.dv <- lisa_s[,1]
shp$lisa_s.p.dv <- p
shp$lisa_s.cl.dv <- as.factor(quadrant)
names(shp)
summary(shp)

'''
# or can merge shapefile with entire data frame from table
nat@data = data.frame(nat@data, nat.lisa[match(nat@data$FIPSNO, nat.lisa$FIPSNO),])
'''

# LISA cluster map, using saddlepoint data

shp@data$id = rownames(shp@data)

#bra@data$values <- rescale(bra$TSLSrho95, to=c(0,1))
shp.points <- fortify(shp, region="id")

'''
# if get error about gpclibPermit not being TRUE, try following:
install.packages("gpclib", type="source")
library(gpclib)
gpclibPermit()
'''

shp.df = join(shp.points, shp@data, by="id")
#shp.df[shp.df==0] <- NA

# Figure 2 (black and white)

plotf <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.cl.dv) + 
  geom_polygon(colour="transparent", fill="white")
plotf

plotf2bw <- plotf + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "black", "grey","",""), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  labs(x="", y="", title="") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plotf2bw

setwd(paste0(path, "/maps", sep=""))
tiff(file="fig2_lisa-bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
print(plotf2bw)
dev.off()

# Figure 2 (color version)

plotf2col <- plotf + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plotf2col

##########################
# other LISA plots; these are not included in manuscript but might be helpful to new users in thinking about their own research

# LISA values

plotv <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.dv) + 
  geom_polygon(colour="transparent", fill="white")
plotv

plotv2 <- plotv + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  labs(x="", y="", title="LISA Values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plotv2


# LISA p-values

plotp <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=lisa_s.p.dv) + 
  geom_polygon(colour="transparent", fill="white")
plotp

plotp2 <- plotp + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, lisa_s.p.dv<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", values=rescale(c(0, .05, .2)), colours=c("darkred", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
plotp2


########################################
#
# CHECKING FOR NON-STATIONARITY 2
# MONTECARLO TEST (RANDOMIZATION)
#
########################################

#library(GWmodel)

# sample code from R manual for GWR (Bivand 2014)
# DM<-gw.dist(dp.locat=coordinates(londonhp))
# bw<-bw.gwr(PURCHASE~FLOORSZ,data=londonhp,dMat=DM, kernel="gaussian")
#See any difference in the next two commands and why?
# res.mont1<-montecarlo.gwr(PURCHASE~PROF+FLOORSZ, data = londonhp,dMat=DM,
# nsim=99, kernel="gaussian", adaptive=FALSE, bw=3000)
# res.mont2<-montecarlo.gwr(PURCHASE~PROF+FLOORSZ, data = londonhp,dMat=DM,
# nsim=99, kernel="gaussian", adaptive=FALSE, bw=300000000000)

coords.so <- coordinates (shp.so)
# create distance matrix based on coords; this can take a while ()
DM.s<-gw.dist(dp.locat=coords.so)
# create bw based on model, distance matrix, and kernel type
# if adaptive=FALSE, finds kernel based on fixed distance band
# if adaptive=TRUE, finds kernel based on number of neighbors (i.e., distance varies)
data.so<-data.frame(shp.so)
names(data.so)
#detach(data.so)
attach(data.so)

bw1<-bw.gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
            data=shp.so,dMat=DM.s, adaptive=TRUE, kernel="bisquare")

# time: bw.gwr command above took less than 5 seconds

# randomisation test (montecarlo) to test null of stationary coefficients
# with large data sets, may want to set nsim low at first, just to check run time
res.mc1.s<-montecarlo.gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
                          data = shp.so,dMat=DM.s, nsim=1000, kernel="bisquare", adaptive=TRUE, bw=bw1)
# also took less than 5 minutes

# check robustness of montecarlo test with higher nsims and different kernel options
# other options for kernel include gaussian, exponenential, tricube, boxcar
# see: help (bw.gwr)
# can also run each of these options with different adaptive settings (TRUE for adaptive or FALSE for fixed)

# results again
res.mc1.s

# write results to files
setwd(paste0(path, "/stationaritytests", sep=""))

write.csv(res.mc1.s, "res.mc1.south.csv") 

setwd(path)
save.image("harbersingram20161123.RData")
#load("gwr-sl20161108.RData")

########################################################################
#####
##### GWR
#####
########################################################################

# CHOOSE KERNEL TYPE AND OPTIMIZE THE BANDWIDTH

# We want to optimize the "neighborhood" for the GWR regressions.
# aic optimization is very slow.  We'll use the cv optimization option.
# Need to choose the kernel shape (bi-square v.s. Gaussian) and whether
# we want fixed or adaptive bandwidth.  We choose here an *adaptive* band
# width and a bisquare kernel: 

#install.packages('spgwr')
#library(spwgr)
#help (gwr.sel)


bw1sl <- gwr.sel(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
                 data=shp.so@data, adapt = TRUE, method = "aic", gweight = gwr.bisquare, coords = coords.so)
# time: gwr.sel command above took about 1 hour

# Notes below are from londonhp example from manual and PSU site
# Optimal adaptive bandwidth for the above case is: approximately 0.32  
# each regression will expand or contract so as to include roughly 32% of sample 

# Similar set up but with a *fixed* bandwidth and Gaussian kernel
# presently the default Gaussian equation

# fixed.bw1 <- gwr.sel(LOGCMEDV ~ CRIM + NOX + LSTAT + DIS, adapt = FALSE,
#	method = "cv", gweight = gwr.Gauss, coords = XYgridtable) 

# Optimal fixed bandwidth: ~ 33698m  This means that the bandwidth for
# each regression using a fixed bandwidth will reach out to include all
# centroids within approximately 34km 

# fixed.bw2 <- gwr.sel(LOGCMEDV ~ CRIM + NOX + LSTAT + DIS, adapt = FALSE,
#	method = "cv", gweight = gwr.bisquare, coords = XYgridtable) 

# Optimal fixed bandwidth: ~ 261213m  This means that the bandwidth for
# each regression using a fixed bandwidth will reach out to include all
# centroids within approximately 261km 


##### FIRST, USING THE ADAPTIVE BISQUARE KERNEL
##### Note that the input for the adapt command is the value generated
##### in the calculation above for adaptive.bw

# help(gwr)
coords.so <- coordinates(shp.so)
data.so <- shp.so@data
detach(data.so)
attach(data.so)

#ptm <- proc.time()
#gwr1 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
#	adapt = 0.1, gweight = gwr.bisquare, coords = coords, hatmatrix = TRUE)      # here setting adapt=.1 manually based on 5 iterations of gwr.sel
#proc.time() - ptm

# parllelize to make much faster
#library(parallel)
#detectCores()  # gives number of cores
#ncores <- detectCores()-1
#cl <- makeCluster(ncores)

#ptm <- proc.time()
#gwr1 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
#            adapt = 0.1, gweight = gwr.bisquare, coords = coords, hatmatrix = FALSE, cl=cl, fit.points=coords)
#proc.time() - ptm


#gwr1.data <- gwr1$SDF
#names(gwr1.data)

# check correlations among GWR coefficients (Wheeler and Tiefelsdorf)
#cor(as.data.frame(gwr1.data[,3:7]))
#rcorr(as.matrix(as.data.frame(gwr1.data[,3:7])), type="pearson")
#pairs(as.data.frame(gwr1.data[,3:7]))


#setwd(paste0(path, "/output", sep=""))  # created new folder 'output' in home directory first before this command
#dir()

#write.csv(rcorr(as.matrix(as.data.frame(gwr1.data[,3:7])), type="pearson")$r, "gwr1.corr.csv")

# fit tests (testing GWR improvement over OLS)
#BFC02.gwr.test(gwr1)
#BFC99.gwr.test(gwr1)
#LMZ.F1GWR.test(gwr1)
#LMZ.F2GWR.test(gwr1)
# other test for varation in local coefficients (stationarity)
#LMZ.F3GWR.test(gwr1)

# graph
#pairs(gwr1.data[,3:7])
#scatterplotMatrix(as.data.frame(gwr1.data[,3:7]))

# graph of scatterplot
#setwd("C:/Users/mi122167/Dropbox/Met Matt/APSA 2016 - Analyzing Diffusion/replication materials/figures")
#tiff(file="figD_gwr_corrs.tif", height=6, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(scatterplotMatrix(as.data.frame(gwr1.data[,3:7])))
#dev.off()


# simple plots
#spplot(gwr1$SDF, "RD90", at=quantile(gwr1$SDF$RD90), col.regions=grey.colors(20), main="Local Coef. for RD90")

# nicer plots
# gen t-values
#gwr1$SDF$RD90_t <- gwr1$SDF$RD90/gwr1$SDF$RD90_se
#gwr1$SDF$PS90_t <- gwr1$SDF$PS90/gwr1$SDF$PS90_se
#gwr1$SDF$MA90_t <- gwr1$SDF$MA90/gwr1$SDF$MA90_se
#gwr1$SDF$DV90_t <- gwr1$SDF$DV90/gwr1$SDF$DV90_se
#gwr1$SDF$UE90_t <- gwr1$SDF$UE90/gwr1$SDF$UE90_se
# try subsetting ggplot; if that does not work; then gen vars at 95% sig
# generate new var with only significant local coefficient (95% level)

#gwr1$SDF$RD90_95 <- gwr1$SDF$RD90
#gwr1$SDF$RD90_95[abs(gwr1$SDF$RD90_t)<1.96] <- 0

#gwr1$SDF$PS90_95 <- gwr1$SDF$PS90
#gwr1$SDF$PS90_95[abs(gwr1$SDF$PS90_t)<1.96] <- 0

#gwr1$SDF$DV90_95 <- gwr1$SDF$DV90
#gwr1$SDF$DV90_95[abs(gwr1$SDF$DV90_t)<1.96] <- 0

#gwr1$SDF$UE90_95 <- gwr1$SDF$UE90
#gwr1$SDF$UE90_95[abs(gwr1$SDF$UE90_t)<1.96] <- 0

#'''
# could also gen new var with only sig loc b at 90% level
#gwr1$SDF$NOM90 <- gwr1$SDF$NOM
#gwr1$SDF$NOM90[abs(gwr1$SDF$NOM_t)<=1.65] <- 0
#'''

# add sig coefficients to shp

#shp$RD90_b <- gwr1$SDF$RD90
#shp$RD90_t <- gwr1$SDF$RD90_t
#shp$RD90_95 <- gwr1$SDF$RD90_95

#shp$PS90_b <- gwr1$SDF$PS90
#shp$PS90_t <- gwr1$SDF$PS90_t
#shp$PS90_95 <- gwr1$SDF$PS90_95

#shp$DV90_b <- gwr1$SDF$DV90
#shp$DV90_t <- gwr1$SDF$DV90_t
#shp$DV90_95 <- gwr1$SDF$DV90_95

#shp$UE90_b <- gwr1$SDF$UE90
#shp$UE90_t <- gwr1$SDF$UE90_t
#shp$UE90_95 <- gwr1$SDF$UE90_95

# set up for ggplot
#shp@data$id = rownames(shp@data)
#bra@data$values <- rescale(bra$TSLSrho95, to=c(0,1))
#shp.points <- fortify(shp, region="id")
#'''
# if get error about gpclibPermit not being TRUE, try following:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()
#'''
#shp.df = join(shp.points, shp@data, by="id")
#shp.df[shp.df==0] <- NA

# try subsetting plot (if this works, then no need to generate RD90_95 or any other variables with only stat sig loc betas)
#plot8 <- ggplot(shp.df) + 
#  aes(long,lat,group=group,fill=RD90_b) + 
#  geom_polygon(colour="transparent", fill="white")
#plot8

#plot8a <- plot8 + 
#  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, abs(RD90_t)>=1.96)) +     # make sure data are attached for this to work
#             colour = 'white', fill = 'black', alpha = .4, size = .3) +
#scale_fill_distiller(name="LISA values", palette = "Reds", 
#                    trans="reverse", 
# breaks = pretty_breaks(n = 5), 
#                   na.value="white", guide="colourbar") + 
#scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
#scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
#labs(x="", y="", title="Local Coef RD90") +
#theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
#axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
#plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
#panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
#panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
#panel.margin=unit(c(0,0,0,0), "lines"),
#plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
#coord_equal(ratio=1) +
#geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
#geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot8a

# send plot to file
# graph of scatterplot
#setwd(paste0(filepath, "/figures", sep=""))
#tiff(file="gwr-rd90.tif", height=6, width=6, units="in", res=600) # could also do pdf, jpeg, bmp, tiff
#print(plot8a)
#dev.off()



# compare to map with RD90_95 to ensure correct
# NOTE: Graphs that do not use subset command are commented out below because subsetting works and is more efficient


####################
# GWR-SL
####################

# review intuition (see Anselin and Rey 2014, pp. 159-162)

# two versions: (A) with crude addition of simple lagged y, and (B) 2SLS approach;
# 2SLS is technically more appropriate, but let's just check and compare how sensitive the results are to the different model specifications

#############################
# (A) crude addition of SL (brute force approach)

# steps
# (1) generate spatial lag of y
# (2) specify gwr-sl model (gwrsl1) by building basic gwr1 and adding lag of y
#############################

# set w
w <- w10s

# step 1: generate spatial lag (Wy); use lag.listw (w needs to be a listw object)
shp.so$HR90lag <- lag.listw(w, HR90)

data.so <- shp.so@data
#detach()
#attach(data)

# step 2: add Wy to gwr to specify gwr-sl
gwrsl1 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + HR90lag, 
              data=shp.so@data, adapt = bw1sl, gweight = gwr.bisquare, coords = coords.so, hatmatrix = TRUE, predictions=TRUE)

names(gwrsl1$SDF)
summary(gwrsl1$SDF$HR90lag)
summary(gwrsl1$SDF$HR90lag_se)

# gen t-values
gwrsl1$SDF$HR90lag_t <- gwrsl1$SDF$HR90lag/gwrsl1$SDF$HR90lag_se
# try subsetting ggplot; if that does not work; then gen vars at 95% sig
# generate new var with only significant local coefficient (95% level)


# add sig coefficients to shp

shp.so$HR90lag_b <- gwrsl1$SDF$HR90lag
shp.so$HR90lag_t <- gwrsl1$SDF$HR90lag_t
#shp$HR90lag_95 <- gwrsl1$SDF$HR90lag_95

# set up for ggplot
shp.so@data$id = rownames(shp.so@data)
#bra@data$values <- rescale(bra$TSLSrho95, to=c(0,1))
shp.so.points <- fortify(shp.so, region="id")
#'''
# if get error about gpclibPermit not being TRUE, try following:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()
#library(maptools)
#'''
shp.so.df = join(shp.so.points, shp.so@data, by="id")
shp.so.df[shp.so.df==0] <- NA

# try subsetting plot (if this works, then no need to generate RD90_95 or any other variables with only stat sig loc betas)
plot1 <- ggplot(shp.so.df) + 
  aes(long,lat,group=group,fill=HR90lag_b) + 
  geom_polygon(colour="transparent", fill="white")
plot1

library(scales)

plot1a <- plot1 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(HR90lag_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  scale_fill_distiller(name="localrho", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  labs(x="", y="", title="Local rho (simple Wy)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot1a

setwd(paste0(path, '/maps', sep=""))

png(file='gwrsl_plain.png', height=6, width=6, units="in", res=300)
print(plot1a)
dev.off()

setwd(path)

save.image("harbersingram20161123.RData")

#######################
# (B) 2SLS 
#######################

# 2 stages:
# (1) in first stage, estimate spatial lag of y (Wy)
# reasonable instruments for Wy are the spatial lags of all Xs
# (review intuition)
# so, this step involves:
# (a) generating spatial lags of Xs
# (b) estimating first stage
# (2) estimate gwr-sl with predicted values from first stage

# stage 1(a): generate lags of Xs, i.e., WX (using lag.listw)

w <- w10s
shp.so$w1RD90 <- lag.listw(w, shp.so$RD90)
shp.so$w1PS90 <- lag.listw(w, shp.so$PS90)
shp.so$w1MA90 <- lag.listw(w, shp.so$MA90)
shp.so$w1DV90 <- lag.listw(w, shp.so$DV90)
shp.so$w1UE90 <- lag.listw(w, shp.so$UE90)

w <- w10.2s
shp.so$w2RD90 <- lag.listw(w, shp.so$RD90)
shp.so$w2PS90 <- lag.listw(w, shp.so$PS90)
shp.so$w2MA90 <- lag.listw(w, shp.so$MA90)
shp.so$w2DV90 <- lag.listw(w, shp.so$DV90)
shp.so$w2UE90 <- lag.listw(w, shp.so$UE90)

data.so <- shp.so@data
detach(data.so)
attach(data.so)

# stage 1(b.1)
# Wy = WX + W2X + W3X ...
HR90ivlag1 <- lm(HR90 ~ w1RD90 + w1PS90 + w1MA90 + w1DV90 + w1UE90)
Wyhat.1 <- HR90ivlag1$fitted.values
# simple Wy
Wy <- lag.listw(w10s, HR90)
plot(Wy, Wyhat.1)

shp.so@data$Wyhat.1 <- Wyhat.1

# stage 1(b.2)
HR90ivlag2 <- lm(HR90 ~ w1RD90 + w1PS90 + w1MA90 + w1DV90 + w1UE90 + w2RD90 + w2PS90 + w2MA90 + w2DV90 + w2UE90)
Wyhat.2 <- HR90ivlag2$fitted.values
plot(Wy, Wyhat.2)

shp.so@data$Wyhat.2 <- Wyhat.2

data.so <- shp.so@data
#detach()
#attach(data)
# alternate stage 1(b) 
#library(parallel)
#detectCores()  # gives number of cores
#ncores <- detectCores()-1
#cl <- makeCluster(ncores)

#ptm <- proc.time()
#ivlag.gwr1 <- gwr(HR90 ~ w1RD90 + w1PS90 + w1MA90 + w1DV90 + w1UE90,
#                      adapt = 0.1, gweight = gwr.bisquare, coords = coords, hatmatrix = FALSE, cl=cl, fit.points=coords)
#proc.time() - ptm

#stopCluster()


# stage 2


gwrsl2.1 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + Wyhat.1,
                data=data.so, adapt = bw1sl, gweight = gwr.bisquare, coords = coords.so, hatmatrix = TRUE)

names(gwrsl2.1$SDF)
summary(gwrsl2.1$SDF$Wyhat.1)
summary(gwrsl2.1$SDF$Wyhat.1_se)

# gen t-values
gwrsl2.1$SDF$Wyhat.1_t <- gwrsl2.1$SDF$Wyhat.1/gwrsl2.1$SDF$Wyhat.1_se
# try subsetting ggplot; if that does not work; then gen vars at 95% sig
# generate new var with only significant local coefficient (95% level)


# add sig coefficients to shp

shp.so$Wyhat.1_b <- gwrsl2.1$SDF$Wyhat.1
shp.so$Wyhat.1_t <- gwrsl2.1$SDF$Wyhat.1_t
#shp$HR90lag_95 <- gwrsl1$SDF$HR90lag_95

# set up for ggplot
shp.so@data$id = rownames(shp.so@data)
#bra@data$values <- rescale(bra$TSLSrho95, to=c(0,1))
shp.so.points <- fortify(shp.so, region="id")
#'''
# if get error about gpclibPermit not being TRUE, try following:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()
#library(maptools)
#'''
shp.so.df = join(shp.so.points, shp.so@data, by="id")
shp.so.df[shp.so.df==0] <- NA

setwd(path)
save.image("harbersingram20161123.RData")

# try subsetting plot (if this works, then no need to generate RD90_95 or any other variables with only stat sig loc betas)
plot2 <- ggplot(shp.so.df) + 
  aes(long,lat,group=group,fill=Wyhat.1_b) + 
  geom_polygon(colour="transparent", fill="white")
#plot2

library(scales)
plot2a <- plot2 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(Wyhat.1_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  labs(x="", y="", title="Local rho (S2SLS, WX)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot2a

setwd(paste0(path, '/maps', sep=""))

png(file='gwrsl_S2SLS-WX.png', height=6, width=6, units="in", res=300)
print(plot2a)
dev.off()

setwd(path)
save.image("harbersingram20161123.RData")



# stage 2, with 2 orders of X


gwrsl2.2 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + Wyhat.2,
                data=shp.so@data, adapt = bw1sl, gweight = gwr.bisquare, coords = coords.so, hatmatrix = TRUE)
# note: in paper, we originally used bw=561, or 0.182
# here, adaptive.bw1 is 0.15, so results may differ slightly from published paper
# still, strategy for case selection remains the same
# as we note in paper, scholars should do several auxiliary checks with different bandwidths and kernel types to check 
# robustness of indicated case selection

names(gwrsl2.2$SDF)
summary(gwrsl2.2$SDF$Wyhat.2)
summary(gwrsl2.2$SDF$Wyhat.2_se)

# gen t-values
gwrsl2.2$SDF$Wyhat.2_t <- gwrsl2.2$SDF$Wyhat.2/gwrsl2.2$SDF$Wyhat.2_se
# try subsetting ggplot; if that does not work; then gen vars at 95% sig
# generate new var with only significant local coefficient (95% level)


# add sig coefficients to shp

shp.so$Wyhat.2_b <- gwrsl2.2$SDF$Wyhat.2
shp.so$Wyhat.2_t <- gwrsl2.2$SDF$Wyhat.2_t
#shp$HR90lag_95 <- gwrsl1$SDF$HR90lag_95

data.so <- shp.so@data
#detach()
#attach(data)
# set up for ggplot
shp.so@data$id = rownames(shp.so@data)
#bra@data$values <- rescale(bra$TSLSrho95, to=c(0,1))
shp.so.points <- fortify(shp.so, region="id")
#'''
# if get error about gpclibPermit not being TRUE, try following:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()
#library(maptools)
#'''
shp.so.df = join(shp.so.points, shp.so@data, by="id")
shp.so.df[shp.so.df==0] <- NA

# try subsetting plot (if this works, then no need to generate RD90_95 or any other variables with only stat sig loc betas)
plot3 <- ggplot(shp.so.df) + 
  aes(long,lat,group=group,fill=Wyhat.2_b) + 
  geom_polygon(colour="transparent", fill="white")
plot3

library(scales)
plot3a <- plot3 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(Wyhat.2_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="localrho", low="lightgrey", high="red", guide="colourbar", na.value="white") +
  scale_fill_gradient(name="localrho", low="lightgrey", high="black", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local rho (S2SLS, W2X)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot3a


plot3b <- plot3 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(Wyhat.2_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="Local b", values=rescale(c(-.8, 0, 1)), colours=c("blue", "white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="localrho", low="lightgrey", high="red", guide="colourbar", na.value="white") +
  #scale_fill_gradient(name="localrho", low="lightgrey", high="black", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local rho (S2SLS, W2X)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 

# break up positive and negative rho coefficients into 2 plots
plot3c <- plot3 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(Wyhat.2_t)>=1.96 & Wyhat.2_b>0)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="localrho", low="lightgrey", high="red", guide="colourbar", na.value="white") +
  scale_fill_gradient(name="localrho", low="lightgrey", high="black", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local rho > 0 (S2SLS, W2X)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot3c

plot3d <- plot3 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.so.df, abs(Wyhat.2_t)>=1.96 & Wyhat.2_b<0)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local b", values=rescale(c(0, 2, 6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  #scale_fill_distiller(name="localrho", palette = "RdYlBu", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="localrho", low="lightgrey", high="red", guide="colourbar", na.value="white") +
  scale_fill_gradient(name="localrho", low="black", high="lightgrey", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local rho < 0 (S2SLS, W2X)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.so, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
#plot3d

setwd(paste0(path, '/maps', sep=""))

png(file='gwrsl_S2SLS-W2X_bw.png', height=6, width=6, units="in", res=300)
print(plot3a)
dev.off()

png(file='gwrsl_S2SLS-W2X.png', height=6, width=6, units="in", res=300)
print(plot3b)
dev.off()

png(file='gwrsl_S2SLS-W2X_pos_bw.png', height=6, width=9, units="in", res=300)
#tiff(file='gwrsl_S2SLS-W2X_pos_bw.tif', height=6, width=9, units="in", res=300)
print(plot3c)
dev.off()

png(file='gwrsl_S2SLS-W2X_neg_bw.png', height=6, width=9, units="in", res=300)
#tiff(file='gwrsl_S2SLS-W2X_neg_bw.tif', height=6, width=9, units="in", res=300)
print(plot3d)
dev.off()

setwd(path)
save.image("harbersingram20161123.RData")

# sorted by local rho from last model, S2SLS, W2X
data.sort <- shp.so@data[order(-shp.so$Wyhat.2_b),]
data.sort <- subset(data.sort, abs(data.sort$Wyhat.2_t)>1.96)
data.sort[1:140,c("NAME","STATE_NAME","FIPSNO","HR90","Wyhat.2_b","Wyhat.2_t")]

temp1 <- data.sort[1:140,c("NAME","STATE_NAME","FIPSNO","HR90","Wyhat.2_b","Wyhat.2_t")]

setwd(paste0(path, "/output", sep=""))
write.csv(temp1, "gwr-sl.top10percent.csv") 

setwd(path)

# compare results across gwrsl1 and both gwrsl2.1 and gwrsl2.2


#####################
#
# GWR-SE (aka, GWR-sEA)
#
#####################

# path <- '/network/rit/lab/ingramlab/geonested'
# or, if running from local computer (not secure shell)
# path <- setwd('A:/geonested')    # where 'A' is mapped network drive
# or, if using directory on disk drive
#path <- 'C:/Users/MI122167/Dropbox/Met Matt/Mixed Methods Paper/PA Submission/replication'

setwd(path)
load("harbersingram20161123.RData")

# shift to non-south data
setwd(paste0(path, '/shapefiles-data', sep=""))

shpns <- readShapeSpatial("nat_nonsouth2.SHP", ID="FIPSNO")
#shp <- readShapeSpatial("NAT.SHP", ID="FIPSNO")
# note: sometimes dbf file will not load if extension is .DBF; if this happens, try renaming to lower case .dbf
#shp.st.so <- readShapeSpatial("tl_2014_us_state_contig_south.shp")
shp.st.ns <- readShapeSpatial("tl_2014_us_state_contig_nonsouth.shp")

###############
#
# Spatial Weights
#
###############

# Here, might only calculate 10 nearest neighbors (10nnb) since that is what is used by Baller et al
# However, will calculate that W (w10), as well as first-order rook (r1) and queen (q1) contiguity matrices 
# row.names(shp)
coords.ns <- coordinates(shpns)
IDs<-row.names(as(shpns, "data.frame"))
w_10nnb<-knn2nb(knearneigh(coords.ns, k=10), row.names=IDs)

# make symmetric
wsym_10nnb <- make.sym.nb(w_10nnb)
wsym_10nnb

# second-order 10nnb w
wsym_10nnb2 <- nblag(wsym_10nnb, maxlag=2)
wsym_10nnb2
# not that this is a list object containing two nb objects
wsym_10nnb2[[1]]  # this is first order w, equivalent to wsym_10nnb
wsym_10nnb2[[2]]  # this is second order w

# if you wanted cumulative lags for all orders within nblag object (here, wsym_10nnb2):
# w.cumul <- nblag_cumul(wsym_10nnb2)

# and make these list objects
listwsym_10nnb <- nb2listw(wsym_10nnb, style="W")
listwsym_10nnb
listwsym_10nnb2 <- nb2listw(wsym_10nnb2[[2]], style="W")
listwsym_10nnb2

# create shortcut to symmetric list objects
w10ns <- listwsym_10nnb
w10ns 

w10.2ns <- listwsym_10nnb2
w10.2ns

# plot
#par(mfrow=c(2,1)) 
#plot(w10, coords, col="red")
#plot(w10.2, coords, col="blue")
#par(mfrow=c(1,1))

#plot(shp)
#plot(w10, coords, col="red", add=T)
#plot(w10.2, coords, col="blue", add=T)

# rook 1
w_nbr1<-poly2nb(shpns, queen=FALSE)
# make symmetric
wsym_nbr1 <- make.sym.nb(w_nbr1)
wsym_nbr1

# second-order 10nnb w
wsym_nbr1.2 <- nblag(wsym_nbr1, maxlag=2)
wsym_nbr1.2

# and make these list objects
listwsym_nbr1 <- nb2listw(wsym_nbr1, style="W")
listwsym_nbr1
listwsym_nbr1.2 <- nb2listw(wsym_nbr1.2[[2]], style="W")
listwsym_nbr1.2
# create shortcuts
wr1ns <- listwsym_nbr1
wr1.2ns <- listwsym_nbr1.2

wr1ns
wr1.2ns

# plot
#par(mfrow=c(2,1)) 
#plot(wr1, coords, col="red")
#plot(wr1.2, coords, col="blue")
#par(mfrow=c(1,1))

# queen 1
w_nbq1<-poly2nb(shpns, queen=TRUE)
# make symmetric
wsym_nbq1 <- make.sym.nb(w_nbq1)
wsym_nbq1

# second-order 10nnb w
wsym_nbq1.2 <- nblag(wsym_nbq1, maxlag=2)
wsym_nbq1.2

# and make these list objects
listwsym_nbq1 <- nb2listw(wsym_nbq1, style="W")
listwsym_nbq1
listwsym_nbq1.2 <- nb2listw(wsym_nbq1.2[[2]], style="W")
listwsym_nbq1.2
# create shortcuts
wq1ns <- listwsym_nbq1
wq1.2ns <- listwsym_nbq1.2

wq1ns
wq1.2ns



# set w
w <- w10ns
#w <- w10.2ns

# step 1: generate spatial lag of residual from OSL model (Wu); use lag.listw (w needs to be a listw object)
detach(data.so)
data.ns <- shpns@data
attach(data.ns)
ols <- lm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90, data = data.ns)
summary(ols)
AIC(ols)
#str(ols)
# will need residual from ols2; ols2$residuals
shpns$u <- ols$residuals
shpns$Wu <- lag.listw(w, shpns$u)

data.ns <- shpns@data
detach(data.ns)
#attach(data)
save.image("harbersingram20161123.RData")

# select bandwidth
coords.ns <- coordinates (shpns)
# create distance matrix based on coords; this can take a while ()
DM.ns<-gw.dist(dp.locat=coords.ns)

# check optimal bw
#bw1<-bw.gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
#            data=shpns,dMat=DM, adaptive=TRUE, kernel="bisquare")

# randomisation test (montecarlo) to test null of stationary coefficients
# with large data sets, may want to set nsim low at first, just to check run time
#res.mc1<-montecarlo.gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90,
#                        data = shpns,dMat=DM, nsim=1000, kernel="bisquare", adaptive=TRUE, bw=bw1)
# also took less than 5 minutes

# check robustness of montecarlo test with higher nsims and different kernel options
# other options for kernel include gaussian, exponenential, tricube, boxcar
# see: help (bw.gwr)
# can also run each of these options with different adaptive settings (TRUE for adaptive or FALSE for fixed)

# results again
#res.mc1

# write results to files
#setwd(paste0(path, "/stationaritytests", sep=""))

#write.csv(res.mc1, "res.mc1.csv") 

setwd(path)


# Note that the general model here is:
# y = XB + u
# where u = Wu + e

# Start with brute force version, adding simple Wu to model and estimate GWR

#select bw
bw1se <- gwr.sel(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + Wu,
                 data=shpns@data, adapt = TRUE, method = "aic", gweight = gwr.bisquare, coords = coords.ns)

gwrse1 <- gwr(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + Wu, 
              data=shpns@data, adapt = bw1se, gweight = gwr.bisquare, coords = coords.ns, hatmatrix = TRUE, predictions=TRUE)

names(gwrse1$SDF)
summary(gwrse1$SDF$Wu)
summary(gwrse1$SDF$Wu_se)

# gen t-values
gwrse1$SDF$Wu_t <- gwrse1$SDF$Wu/gwrse1$SDF$Wu_se
# try subsetting ggplot; if that does not work; then gen vars at 95% sig
# generate new var with only significant local coefficient (95% level)


# add sig coefficients to shp

shpns$Wu_b <- gwrse1$SDF$Wu
shpns$Wu_t <- gwrse1$SDF$Wu_t
#shp$HR90lag_95 <- gwrsl1$SDF$HR90lag_95

# set up for ggplot
shpns@data$id = rownames(shpns@data)
#if need to rescale, i.e., normalize, use:
# shp@data$values <- rescale(shp$[var], to=c(0,1))
shpns.points <- fortify(shpns, region="id")
#'''
# if get error about gpclibPermit not being TRUE, try following:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()
#library(maptools)
#'''
shpns.df = join(shpns.points, shpns@data, by="id")
shpns.df[shpns.df==0] <- NA

# try subsetting plot (if this works, then no need to generate RD90_95 or any other variables with only stat sig loc betas)
plot10 <- ggplot(shpns.df) + 
  aes(long,lat,group=group,fill=abs(Wu_b)) + 
  geom_polygon(colour="transparent", fill="white")
plot10

library(scales)

plot10a <- plot10 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shpns.df, abs(Wu_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="Local lambda", values=rescale(c(0.4, 1, 1.6)), colours=c("white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="locallamda", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  labs(x="", y="", title="Local lambda (absolute value, simple Wu)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.ns, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
plot10a

setwd(paste0(path, '/maps', sep=""))

png(file='gwrse_plain_abs_color.png', height=6, width=6, units="in", res=300)
#tiff(file='gwrse_plain_color.tif', height=6, width=6, units="in", res=300)
print(plot10a)
dev.off()

# b&w version of gwr-se plain
plot10b <- plot10 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shpns.df, abs(Wu_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="Local lambda", values=rescale(c(0.4, 1, 1.6)), colours=c("white", "black"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="locallamda", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  labs(x="", y="", title="Local lambda (absoluate value, simple Wu)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.ns, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
plot10b

setwd(paste0(path, '/maps', sep=""))


png(file='gwrse_plain_abs_bw.png', height=6, width=6, units="in", res=300)
#tiff(file='gwrse_plain_bw.tif', height=6, width=6, units="in", res=300)
print(plot10b)
dev.off()

# set pos and neg values
plot10.5 <- ggplot(shpns.df) + 
  aes(long,lat,group=group,fill=Wu_b) + 
  geom_polygon(colour="transparent", fill="white")
plot10.5

library(scales)

plot10.5a <- plot10.5 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shpns.df, abs(Wu_t)>=1.96)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="Local lambda", values=rescale(c(-1.75, 0, 1)), colours=c("blue", "white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="locallamda", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  labs(x="", y="", title="Local lambda (simple Wu)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.ns, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
plot10.5a

setwd(paste0(path, '/maps', sep=""))

png(file='gwrse_plain_color.png', height=6, width=6, units="in", res=300)
#tiff(file='gwrse_plain_color.tif', height=6, width=6, units="in", res=300)
print(plot10.5a)
dev.off()

# b&w version of pos and neg
plot10.5b <- plot10.5 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shpns.df, abs(Wu_t)>=1.96 & Wu_b>0)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local lambda", values=rescale(c(-1, 0, 1)), colours=c("blue", "white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="locallamda", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  scale_fill_gradient(name="locallambda", low="lightgrey", high="black", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local lambda > 0 (simple Wu)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.ns, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
plot10.5b

plot10.5c <- plot10.5 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shpns.df, abs(Wu_t)>=1.96 & Wu_b<0)) +     # make sure data are attached for this to work
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="Local lambda", values=rescale(c(-1, 0, 1)), colours=c("blue", "white", "red"), na.value="white", guide="colourbar") +
  #scale_fill_distiller(name="locallamda", palette = "Reds", breaks = pretty_breaks(n = 5), trans="reverse", na.value="white", guide="colourbar") + 
  scale_fill_gradient(name="locallambda", low="black", high="lightgrey", guide="colourbar", na.value="white") +
  labs(x="", y="", title="Local lambda < 0 (simple Wu)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1), # make title bold and add space
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.margin=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  #geom_polygon(data=shp, aes(long,lat, group=group), colour="grey50", fill=NA) +
  geom_polygon(data=shp.st.ns, aes(long,lat, group=group), colour="grey10", fill=NA) 
# geom_polygon(data=subset(shp.df, id="42041"), aes(long,lat, group=group), colour="blue")
plot10.5c


setwd(paste0(path, '/maps', sep=""))

#png(file='gwrse_plain_pos_bw.png', height=6, width=9, units="in", res=300)
tiff(file='gwrse_plain_pos_bw.tif', height=6, width=9, units="in", res=300)
print(plot10.5b)
dev.off()

#png(file='gwrse_plain_neg_bw.png', height=6, width=9, units="in", res=300)
tiff(file='gwrse_plain_neg_bw.tif', height=6, width=9, units="in", res=300)
print(plot10.5c)
dev.off()

setwd(path)

save.image("harbersingram20161123.RData")

# sorted by local lambda
summary(shpns$Wu)
# descending order
data.sort <- shpns@data[order(-shpns$Wu_b),]
data.sort <- subset(data.sort, abs(data.sort$Wu_t)>1.96)
data.sort[1:166,c("NAME","STATE_NAME","FIPSNO","HR90","Wu_b","Wu_t")]

temp1 <- data.sort[1:166,c("NAME","STATE_NAME","FIPSNO","HR90","Wu_b","Wu_t")]

setwd(paste0(path, "/output", sep=""))
write.csv(temp1, "gwr-se.top10percent.csv") 

# ascending order
data.sort <- shpns@data[order(shpns$Wu_b),]
data.sort <- subset(data.sort, abs(data.sort$Wu_t)>1.96)
data.sort[1:166,c("NAME","STATE_NAME","FIPSNO","HR90","Wu_b","Wu_t")]

temp1 <- data.sort[1:166,c("NAME","STATE_NAME","FIPSNO","HR90","Wu_b","Wu_t")]

setwd(paste0(path, "/output", sep=""))
write.csv(temp1, "gwr-se.bottom10percent.csv") 

setwd(path)


#end
